/**

This file is part of MaCI/GIMnet.

MaCI/GIMnet is free software: you can redistribute it and/or modify it 
under the terms of the GNU Lesser General Public License as published 
by the Free Software Foundation, either version 3 of the License, or 
(at your option) any later version.

MaCI/GIMnet is distributed in the hope that it will be useful, but WITHOUT 
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public 
License for more details.

You should have received a copy of the GNU Lesser General Public 
License along with GIMnet. (See COPYING.LESSER) If not, see 
<http://www.gnu.org/licenses/>.

**/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int choldcb(double *pR, double *pA, int N, int  p);
void cholslb(double *A, double *b, double *x, int N, int p);
int choldc(double *pR, double *pA, int N);
void cholsl(double *A, double *b, double *x, int N);
void ldumat(double *,double *,int,int);
void ldvmat(double *,double *,int );
int qrbdv(double *,double *,double *,int ,double *,int);
int svduv(double *,double *,double *,int,double *,int);
typedef struct complex_t {
/* Surface plane parameters */
	double	real;
	double	imag;
} complex;

int complex_next_root(double *a,double *a1, complex *v, int *p,int N,int cyc);
void horner(double *a, double *b, double x, int N);
void complex_horner(double *a, double *b, complex x, int N);
double evaluate_poly(double *a, double v, int N);
double complex_evaluate_poly(double *a, complex v, complex *f,int N);
void find_positive_real_roots(double *a,double *a1, double *b, double *x, int N,int cyc,int *k);
void solve_quadratic(complex *x1, complex *x2,double *a);
double power_iter(double *A,double *v, double *x,int N);
void back_substitution(double *A, double *b, double *x, int N);
double find_largest_value(double *v,int N);
int calc_qr_light(double *A,double *R,int N);
void calc_qr(double *Q, double *P, double *R, double *v,double *w,int M,int N);
void find_smallest_eigenvalue(double *A,double *Q,double *R,double *iA, double *v, double *x, int N,int cyc);
